subroutine rgp_fortran_der_hc_aabb( &
       no, & ! (in)    number of occupied orbitals
       nv, & ! (in)    number of unoccupied orbitals
    Tiajb, & ! (in)    two electron integrals in 1122/aabb order
      Via, & ! (in)    new rows
      Zia, & ! (in)    derivative contracted new rows
      Ujb, & ! (in)    new columns
      Yjb, & ! (in)    derivative contracted new columns
      Pab, & ! (in)    pairing matrix of unoccupied sites (indices are alpha beta)
      Qab, & ! (in)    derivative contracted pairing matrix of unoccupied sites
      Ckl, & ! (in)    derivative contracted inverse matrix (indices are beta alpha)
      Clk, & ! (in)    derivative contracted transpose of inverse matrix
     VBia, & ! (in)    contraction of new rows with inverse
     VCia, & ! (in)    part of the derivative of the contraction of new rows with inverse
     ZBia, & ! (in)    part of the derivative of the contraction of new rows with inverse
     BUjb, & ! (in)    contraction of new cols with inverse
     CUjb, & ! (in)    part of the derivative of the contraction of new cols with inverse
     BYjb, & ! (in)    part of the derivative of the contraction of new cols with inverse
    VBUab, & ! (in)    contraction of new rows, new cols, and inverse
     TBab, & ! (in)    contraction of tei and inverse
    TVBjb, & ! (in)    contraction of tei, new rows, and inverse
   dVBUab, & ! (out)   derivative contracted contraction of new rows, new cols, and inverse
     TCab, & ! (out)   contraction of tei and inverse
    TBUia, & ! (out)   contraction of tei, new cols, and inverse
   retval)   ! (inout) running total for the energy

implicit none

integer,      intent(in)    ::     no
integer,      intent(in)    ::     nv
real(kind=8), intent(in)    ::  Tiajb(no,nv,no,nv)
real(kind=8), intent(in)    ::    Via(no,nv) ! second index labels which row
real(kind=8), intent(in)    ::    Zia(no,nv) ! second index labels which row
real(kind=8), intent(in)    ::    Ujb(no,nv) ! second index labels which col
real(kind=8), intent(in)    ::    Yjb(no,nv) ! second index labels which col
real(kind=8), intent(in)    ::    Pab(nv,nv) ! alpha beta
real(kind=8), intent(in)    ::    Qab(nv,nv) ! alpha beta
real(kind=8), intent(in)    ::    Ckl(no,no) ! beta alpha
real(kind=8), intent(in)    ::    Clk(no,no) ! alpha beta
real(kind=8), intent(in)    ::   VBia(no,nv) ! second index labels which row
real(kind=8), intent(in)    ::   VCia(no,nv) ! second index labels which row
real(kind=8), intent(in)    ::   ZBia(no,nv) ! second index labels which row
real(kind=8), intent(in)    ::   BUjb(no,nv) ! second index labels which col
real(kind=8), intent(in)    ::   CUjb(no,nv) ! second index labels which col
real(kind=8), intent(in)    ::   BYjb(no,nv) ! second index labels which col
real(kind=8), intent(in)    ::  VBUab(nv,nv)
real(kind=8), intent(in)    ::   TBab(nv,nv)
real(kind=8), intent(in)    ::  TVBjb(no,nv)
real(kind=8), intent(out)   :: dVBUab(nv,nv)
real(kind=8), intent(out)   ::   TCab(nv,nv)
real(kind=8), intent(out)   ::  TBUia(no,nv)
real(kind=8), intent(inout) :: retval

integer :: a, b, i, j, k, l

include 'formic/fqmc/rgp_fortran_interface.fpp'

! compute dVBUab
dVBUab = 0.00d+00
call rgp_fortran_build_vbuab_aabb(no, nv, Zia, BUjb, dVBUab)
call rgp_fortran_build_vbuab_aabb(no, nv, Via, CUjb, dVBUab)
call rgp_fortran_build_vbuab_aabb(no, nv, Via, BYjb, dVBUab)

! compute TCab
call rgp_fortran_build_tbab_aabb(no, nv, Tiajb, Clk, TCab)

! compute TBUia
call rgp_fortran_build_tbuia(no, nv, Tiajb, BUjb, TBUia)

! term 6 and 16 (derivative hits non-Bji tensors)
do b = 1,nv
do a = 1,nv
  retval = retval + TBab(a,b) * ( Qab(a,b) - dVBUab(a,b) )
enddo
enddo

! term 6 and 16 (derivative hits Bji tensor)
do b = 1,nv
do a = 1,nv
  retval = retval + TCab(a,b) * ( Pab(a,b) - VBUab(a,b) )
enddo
enddo

! term 12
retval = retval + ddot(no*nv, TVBjb, 1, CUjb, 1)
retval = retval + ddot(no*nv, TVBjb, 1, BYjb, 1)
retval = retval + ddot(no*nv, TBUia, 1, ZBia, 1)
retval = retval + ddot(no*nv, TBUia, 1, VCia, 1)

end subroutine rgp_fortran_der_hc_aabb
